HDR Imaging and Display

How can those who live in the light of the day possibly comprehend the depths of the night?
--Nietzsche
Author: Zhen Tong

In this project, we refer to the 2002 paper for HDR image processing. Our contribution includes implementation:

LDR Fusion
Bilateral Filter in numpy
Tone Mapping in YUV space

Step 1. Merge LDR RAW into 32-bit HDR Image

Given a list of RAW image of different light exposure, putting them through the first few process of isp, Dead Pixel Correction, Black Level Compensation, Anti Aliasing Filter, and Auto White Balance and Gain Control. Then we can get on the Raw Exposure Fusion process.

\begin{matrix} I_{i} = f (t_{i} X) \\ I_{l i n_{i}} = f^{- 1} (I_{i}) \\ f (x) = x^{1 / γ} \end{matrix}

$t$ jpg $t_i$ $f$ $x_{ij}\in X_k$ $k$ .

w_{i j} = \exp (- 4 \frac{(x_{i j} - median(X))^{2}}{{median(X)}^{2}})

$w_{ij}$ $t_k$ $(i, j)$ $w_{ij}\rightarrow1$ $t_k$ , and should be fully used in the merging step.

$X$ in HDR, we can minimize the weighted least-squares objective function in log-space:

\begin{matrix} min_{X_{i j}} \sum_{k} w_{i j}^{k} (\log (I_{i, j}^{k}) - \log (t^{k} X_{i j}^{k})^{2}) \\ solve by gradiant=0: {\hat{X}}_{i j} = \exp (\frac{\sum_{k} w_{i j}^{k} (\log (I_{i, j}^{k}) - \log (t^{k}))}{\sum_{k} w_{i j}}) \end{matrix}

Output

get_fusion_weights() $w_{ij}=1$ "bell" shape $w_{ij}=1$ are close to the median value of the distribution).

Output Weight as Image:

merging LDR Exposurs.png $w_{ij}\in[0, 1]$ , we can observe that the the short time exposure, like picture 0, has large weight at the dark place, while the long time exposure, like picture 7, has large weight at the bright place. This is because according to the assumption, the accurate pixel will concentrate at the middle point of the exposure time.

Step 2. Demosaic the Fused Raw Data & Save `.EXR` File

After the merging, using the Malvar Color Filter Array Interpolation and output the rbg image. Using OpenEXR, Imath


x
def save_hdr_image_to_exr(rgb_matrix, file_path):
    # Convert the RGB matrix to a 32-bit float array
    hdr_data = np.array(rgb_matrix, dtype=np.float32)

    # Create an OpenEXR image header with the dimensions of the matrix
    header = OpenEXR.Header(hdr_data.shape[1], hdr_data.shape[0])

    # Set the channels in the header
    header['channels'] = {
        'R': Imath.Channel(Imath.PixelType(Imath.PixelType.FLOAT)),
        'G': Imath.Channel(Imath.PixelType(Imath.PixelType.FLOAT)),
        'B': Imath.Channel(Imath.PixelType(Imath.PixelType.FLOAT))
    }

    # Create an OpenEXR image object
    exr_image = OpenEXR.OutputFile(file_path, header)

    # Convert the RGB matrix to a planar configuration (separate R, G, B channels)
    r_channel = hdr_data[:, :, 0].tobytes()
    g_channel = hdr_data[:, :, 1].tobytes()
    b_channel = hdr_data[:, :, 2].tobytes()

    # Write the channel data to the EXR image
    exr_image.writePixels({'R': r_channel, 'G': g_channel, 'B': b_channel})

    # Close the EXR image file
    exr_image.close()

Output

Fuse seven pictures into one.

Step 3. Tone Mapping with Bilateral Filter

$[0, 255]$ , or else, the RGB image have the overexposure problem if output directively, and the image above is also not good enough the local information lost.

The idea of Tone Mapping is separate the the base and detail of the Y channel in YUV, and using the bilateral filter to remain the detail.

$G$ is a Gaussian function with mean = 0.

\begin{matrix} B F [I]_{p} = \frac{1}{W_{p}} \sum_{q \in S} G_{σ_{s}} (| p - q |) G_{σ_{r}} (I_{p} - I_{q}) I_{q} \\ W_{p} = \sum_{q \in S} G_{σ_{s}} (| p - q |) G_{σ_{r}} (I_{p} - I_{q}) \end{matrix}

Implement BF

Using the vectorization of the numpy computation, we can speed up the computation. We can first get all the indices of the kernel computation locations using np.meshgrid, and then do the kernel computation all at once.


xxxxxxxxxx
h, w = img.shape[0], img.shape[1]
img_pad = np.pad(img, radius+1)
print(img_pad.shape)
ys, xs = np.meshgrid(np.arange(h), np.arange(w))
P = img_pad[ys+radius, xs+radius]
W = np.zeros_like(img.T)
base = np.zeros_like(img.T)
for i in tqdm(range(2*radius-1)):
    for j in tqdm(range(2*radius-1)):
        Q = img_pad[ys+i, xs+j]
        val = P-Q
        gr = 1/(sig1*np.sqrt(2*np.pi)) * np.exp(-0.5*(val**2)/(sig1**2))
        dis = np.sqrt((i-radius)**2+(j-radius)**2 )
        gs = 1/(sig2*np.sqrt(2*np.pi)) * np.exp(-0.5*(dis**2)/(sig2**2))
        grgs = gr*gs 
        base += grgs*Q
        W += grgs
base/=W

Output

According to the paper:

we only use a two scale decomposition, where the “base” image is computed using bilateral filtering, and the detail layer is the division of the input intensity by the base layer.

First we change the output of CFA interpolation directly from RBG into YUV and then use the Bilateral Filter to get the base, and get the detail using:

Y_{d e t a i l} = \frac{Y}{Y_{b a s e}}

Then, we can get the new Y channel with:

\begin{matrix} Y_{n e w} = f_{n o r m a l i z e} (\log (Y_{b a s e}) + Y_{d e t a i l}) \end{matrix}

$f_{normalize}$ $[16,235]$ .

Apply the new Y channel data, and recompose with the UV channel, convert it into RGB, we can get the output.

Then we can go through the rest of the ISP process. (The color of the sky slightly changes)

Compared with Gamma without Tone Mapping:

Using the rest of the ISP pipeline (Edge Enhancement, Brightness/Contrast Control, False Color Suppresion, and Hue/Saturation control), we can get the final output.

output

The green color on the left is because the trees reflect green light. The yellow on the right is because the stone reflects yellow light. So hdr can save these color details after preserving the texture of the leaves

How to Run

You should have the file structure like this:


xxxxxxxxxx
/mnt/DDA4310/src$ tree
.
├── data        
│   ├── data.rar
│   ├── set01
│   │   ├── DSC00049.ARW
│   │   ├── DSC00049.JPG
│   │   ├── DSC00050.ARW
│   │   ├── DSC00050.JPG
│   │   ├── DSC00051.ARW
│   │   ├── DSC00051.JPG
│   │   ├── DSC00052.ARW
│   │   ├── DSC00052.JPG
│   │   ├── DSC00053.ARW
│   │   └── DSC00053.JPG
│   ├── set02
│   │   ├── DSC00113.ARW
│   │   ├── DSC00113.JPG
│   │   ├── DSC00114.ARW
│   │   ├── DSC00114.JPG
│   │   ├── DSC00115.ARW
│   │   └── DSC00115.JPG
│   └── set03
│       ├── DSC00163.ARW
│       ├── DSC00163.JPG
│       ├── DSC00164.ARW
│       ├── DSC00164.JPG
│       ├── DSC00165.ARW
│       ├── DSC00165.JPG
│       ├── DSC00166.ARW
│       ├── DSC00166.JPG
│       ├── DSC00167.ARW
│       ├── DSC00167.JPG
│       ├── DSC00168.ARW
│       ├── DSC00168.JPG
│       ├── DSC00169.ARW
│       └── DSC00169.JPG
└── src
    ├── BayerDomainProcessor.py
    ├── HW03.md
    ├── HW03.pdf
    ├── RGBDomainProcessor.py
    ├── YUVDomainProcessor.py
    └── notebook.ipynb

Environment


xxxxxxxxxx
cv2
numpy
matplotlib
OpenEXR
Imath

Args	Value	Help
--draw-intermediate		Draw the intermediate pictures for the intermediate stage
--set-dir	../data/set03
--output-dir	data	The final output image you want to save
--itmd-dir	itmd	The intermediate path for stage you want to save
--fast	store true	Use the faster version of the Bilateral Filter
--kernel-size	5	The kernel size of the BF

Run as ipynb Notebook

For the status check, we recommend you run the notebook.ipynb for detail output image of each state.

Experiments

See the running log at Runing Log

References

[知乎]亮度响应与HDR基础: https://zhuanlan.zhihu.com/p/23981690
Brown University Project 5: High Dynamic Range: https://cs.brown.edu/courses/cs129/2011/asgn/proj5/
Standford Slides: https://web.stanford.edu/class/ee367/slides/lecture6.pdf
Fr´edo Durand and Julie Dorsey Fast Bilateral Filtering for the Display of High-Dynamic-Range Images. Siggraph 2002: https://stanford.edu/class/ee367/reading/FastBilateralFilteringforHDR.pdf
CMU Lecture Notes: http://graphics.cs.cmu.edu/courses/15-463/2018_fall/lectures/lecture6.pdf
Brown Project Guide: https://cs.brown.edu/courses/csci1290/labs/lab_bilateral/index.html

HDR Imaging and Display

Step 1. Merge LDR RAW into 32-bit HDR Image

Output

Step 2. Demosaic the Fused Raw Data & Save .EXR File

Output

Step 3. Tone Mapping with Bilateral Filter

Implement BF

Output

How to Run

Environment

Run as ipynb Notebook

Experiments

References

Step 2. Demosaic the Fused Raw Data & Save `.EXR` File